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Dynamics of fluctuation of the top location of a sandpile 

Chiyori Urabe* 

Graduate School of Human and Environmental, Kyoto University, Kyoto 606-8501 

We investigate the fluctuation of the top location of a sandpile numerically using the two- 
dimensional discrete elements method. We feed particles to a sandpile at a fixed time interval 
\ and calculate power spectra from the time series of the top location. We find that the power 

spectra are approximated as power functions, and their exponents increase to —1 through a 
plateau as the time interval decreases. For small time interval, avalanches occur continuously 
either on the left or right side slope of a sandpile, and the slope on which avalanches take 
place switches intermittently. The long time fluctuation of the top location corresponds to 
the switchings. For the time series of the switchings, we discuss the relation between the 
' power spectrum and the distribution of waiting times based on analytic theories. 
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1. Introduction 



Dense systems of granular materials exhibit solid-like and fluid-like behaviors, 1 4 and many 



researches are devoted on them. For the static phase, the spatial distribution of stress on the 
bottom of a sandpile was measured in experiments, and it is known that its functional form 
depends on the history of the formation of the sandpile. 5-8 For the fluid phase, steady flows 
and avalanches are investigated from several perspectives. 9-18 Fluidized phase appears in a 
localized layer near the surface, and the particle velocity under the layer obeys an exponential 
function of the depth from the surface. 19 It is known for granular flows in a pipe that the 
temporal power spectra of particle density obey power laws. 20-31 

In systems of a sandpile, feeding particles at small feed rate, the surface of a sandpile is 
kept in solid states except when avalanches occur intermittently, and continuous flows appear 
as the feed rate increases. However, even in the case that the feed rate is rather large, it is 
infrequent that the whole surface of a sandpile is kept in the fluid phase, and the states of 
the surface change temporally or spatially. 32 The time evolution of a sandpile is caused by 
complicated interactions between avalanches and the shape of the surface. 

A sandpile typically becomes mountain shaped with a top as it grows. 33,34 Because parti- 
cles on the surface run down from the top to the foot of a sandpile, the top plays the role of a 
singular point in an average flow of particles. In the formation process of a sandpile, the top 
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location moves with time and determines the global flow on the surface. In this paper, we focus 
on the fluctuation of the top location to characterize the long time evolution of a sandpile. 
We carry out two-dimensional numerical simulations and measure the power spectrum of the 
top location. We find that the power spectrum obeys a power function, and that the exponent 
depends on the feeding rate. 

The organization of the remainder of this paper is as follows. In §2, we explain the method 
of our numerical simulations. In §3, we investigate the power spectrum of the top location 
and relations between its exponent and the feed rate of particles. In §4, we discuss the origin 
of this power law. Finally, we draw conclusions in §5. 

2. Method Of Simulations 

In order to simulate polydispersed circular particles, we adopt a two-dimensional Discrete 
Elements Method (DEM). 35 In DEM, we assume that particles are cohesionless, and the mass 
density of a particle per unit area, po, is constant. The linear spring model is adopted to 
describe the repulsion of two particles in contact, and the viscous forces and Coulomb slip 
are assumed to act on them. We assume that the ith particle has the radius r*, the mass 
rrii = nrfpo and the momentum of inertia Ii = rriirf /2. The center of mass and the angular 
velocity of the ith particle represent Xj and u>i, respectively. The equations of motion of the 
ith particle are given by 



where njj and tij represent the normal and the tangential unit vectors at the contact point 
between the jth particle and the ith particle, and g is the acceleration of gravity, as depicted 




(1) 
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Fig. 1. The ith particle and the jth particle are in contact. 
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Fig. 2. Setup of a sandpile in our simulations. x top represents the horizontal top location. Fed particles 
are released at the height H from the top. 



in Fig. 1. The normal contact force Fn is denned by 

F$ = F$e{-F$), (2) 

and 

F% = -k n m l r j (ri + rj - |xj - Xj |) - rjm^riij ■ (±i - , (3) 

where Q(x) is the Heaviside function, and m % r is the reduced mass ml- = J™*™! . . ©(— x) is 
introduced so that the contact force is repulsive. The tangential contact force F\ 3 is defined 
by 

F l t j = k t r4 j ui j . (4) 
u t J is obtained by integrating the equation 

4 j = -(( Xi - Xj ) -tij + rm + W )e(^| - |^'|), (5) 

where u t J is zero when the particles are not in contact (i. e. for |xj — Xj| > r j + ri). Here, 
Coulomb frictional coefficient fi is assumed to be 0.5. We express the maximum diameter and 
the maximum weight of a particle as d and m. The distribution function of the diameters 
is uniform in the range between 0.8d and d. We assume that the spring constants are k n = 
1.0 x 10 4 mg/d in the normal direction, and kt = 2.0 x W 3 mg/d in the tangential direction, and 
that the viscosity is rj = 1.0 x 10 2 y/g/d. The coefficient of restitution in our model is about 
0.2 for a head-on particles collision. We adopt the second-order Adams-Bashforth method for 
time-integration with the time interval At = 1.0 x 10~ 3 \fdfg. 

We assume that a particle is in a sandpile if the particle contacts other particles, and 
the top location of the sandpile is defined as the center of mass of the highest particle in the 
sandpile (Fig. 2). The floor under a sandpile is a horizontal array of 80 fixed particles with 
diameter d. We introduce the x coordinate along the floor and define the origin at the center 
of the floor. 

To investigate the fluctuation of the top in the formation process of a sandpile, we drop 
particles with the time interval T to it. The particles are released at a position just above 
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the center of the floor and its height is H from the top location as shown in Fig. 2. We first 
make a sandpile grow until it covers the floor and use this sandpile as an initial state. Because 
particles run off the edges of the floor with finite length, the size of a sandpile is maintained 
almost constant. 

After the time series of the top location xt op (t) reaches a statistical stationary state, 
we calculate its power spectrum S(f) with respect to frequency /. To calculate the power 
spectrum S(f) from the time series, we divide it into M time series with a time interval T^ s \ 
the rath power spectrum S m (fj) is defined by 

^(/,)^lX>Sn^'| 2 (6) 

71=1 

where x£J n = x top {(m + f )T«) and fj = ^. 

We introduce the power spectrum <!>(/) as the average of the M power spectra, 

M-\ 

s ^ E iE S M)- (7) 

m=0 

3. Results 

We measure xt op for various values of the time interval T and the height H. We change 
H in the range 20d < H < HOd. If H is sufficiently large beyond this range, the impact of a 
dropped particle is large and collapses the top shape of a sandpile into a caldera. 36 Figures 
3(a) and 3(b) shows the time series of xt op (t) obtained from simulations with T = 2\J d/g (a) 
and T = 80 \fdfg~ (b). The top fluctuates frequently in the case of small T, on the other hand, 
in the case of large T , the top almost stays for long time in comparison with T because the 
motion of particles induced by the impact of a fed particle ceases before the next particle is 
dropped. Figure 4 is the power spectra of the time series, S(f), which is calculated using eqs. 
(6) and (7) with = N = 10000 and M = 10. We find that S(f) behaves as a power-law, 
and its exponent depends on the value of T. At T = 2^/d/g, S( f) approximately obeys 1/ / 
law. At T = 80^ d/g, the exponents of S(f) are smaller than —1. 

We investigate the dependence of the exponent of S(f), a, on T and H. For the range 
5/yW < f < 1/(2T), we calculate a from the double logarithmic plot of S(f) in the least 
square method, a is insensitive to H as shown in Fig. 5(a). In contrast, Fig. 5(b) indicates 
that a strongly depends on T for small T and approaches — 1 as T decreases. As a approaches 
— 1, the power spectrum is approximated as a power function with a high degree of accuracy 
as indicated with the error bars. In the range 10 y/ d/g < T < 60^ d/g, a is approximately 
a constant —1.43 ± 0.03, although a decreases as T increases beyond 60\/d/g. In the case of 
large T, the error bars in Figs. 5 are large because the range of frequency used for fitting is 
small. In the region of higher frequency than 1/T, the power function with the same exponent 
is not best fit with S(f). If we fit S(f) with a power function in this range of high frequency, 
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Fig. 3. Time series of the top location x top (t)- We use H = 110c?, T = I^Jdjg (a) and T — SQy'd/g 
(b). The vertical axes are time t. We note that the same number of particles are fed in (a) and 
(b). 



its exponent changes from that indicated in Figs. 5, as shown with the data of T = 80a/ d/g 
in Fig. 4. 

We mainly focus on the case of small T because we are, in particular, interested in the case 
that a is close to —1. The displacement of the top location xtop is caused by avalanches, and 
avalanches occur on either slope at almost all times in this case. Although the instantaneous 
magnitude of avalanches is characterized by the kinetic energy of the particles, S(f) and the 
power spectrum of the kinetic energy differ in functional form as mentioned below. Eliminating 
the narrow region — d < x < d in the center of a sandpile, we divide the sandpile into the left 
part and the right part with respect to x = 0. We measure the kinetic energies of the left and 
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right parts, Ki(t) and K r (t), respectively. Using the same definition in eqs. (6) and (7), the 
power spectra of the time series of Ki{t) and K r {t) are calculated with = 1.0 x ltfy/d/g, 
N = 1.0 x 10 5 and M = 10. For T = and H = 20d, the power spectra of both K t (t) 

and xto P (t) are shown in Fig. 6. The power spectrum of K r (t) is similar to that of Ki(t). 
Because the power spectra of Ki(t) and K r {t) are Lorentzian-like, avalanches seem to occur 
at random. 

From the results of numerical simulations, it is found to be rare that avalanches occur 
simultaneously on the left and right slopes of a sandpile. We refer the states that avalanches 
occur on the left and right slopes as left mode and right mode, respectively. To investigate 
switchings between the both mode, we define the binarized time series, 

+1, for Ki(t) < K r (t), 



K(t) = I 



(8) 

-1, for Ki(t) > K r (t). 



The sign of K(t) represents the side on which avalanches occur mainly at time t. The switchings 
are well defined by K(t) in the case that T is small. However, as T increases, it is difficult 
to define the switchings because avalanches occur at intervals, and the time intervals between 
avalanches are comparable to the time scale of switchings. We find that the time series of 
K(t) is similar to xtop{t) for small T. The power spectrum of the time series K(t) is shown 
in Fig. 7 for T = 2y 'd/g and H = 20d. The power spectrum of K(t) is approximated as a 
power law with the exponent of —1 in the long time scale. The exponent is approximately the 
same as that of the top location. Investigating the conditional probability of K(t) = —1 for a 
given xt p, this probability increases with the value of xt op - Therefore, in each mode, the top 
location xt op {t) is mainly in the opposite side on which avalanches occur. Thus the fluctuation 
of xtop corresponds to the switching between the two modes, but not to the fluctuation of the 
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Fig. 5. (a) Dependence of a on If and (b) dependence of a on T. The error bars correspond to the 
standard deviations. 



magnitude of avalanches. 
4. Discussion 

For the binarized time series such as K(t) defined by eq. (8), it is known from an analytical 
theory that its spectrum is expressed as a power function if the waiting time has a power-law 
distribution and each interval is independent. 37,38 Here, the waiting time r is defined as a 
time interval between neighboring switchings in the binarized time series. We assume that the 
probability density of r, p(r), is the abrupt-cutoff power law, 



p(r) = < 



ct D , for a < t < b, 

(9) 

0, otherwise, 
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Fig. 6. Power spectra of the time series of Xt op (t) and K[(t) for T = 2-v/ d/g and H = 20d. The line 
with a slope is drawn for reference. 
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Fig. 7. Power spectra of the time series of x top (t) and K(t) for T = 2^ d/g and H = 20d. 

where the constants a and b are sufficiently small and large respectively, and c is the normal- 
ization constant. In the range of 1/6 <C / 1/a, the power spectrum of this binarized time 
series, Sb(f), is given approximately by 

s b (f) ~ /- (D+3) (10) 

for —3 < D < —1. In the case of small T, p(r) is expected to be a power function with D = —2 
because the exponent of the power spectrum of K(t) is approximately —1 in our simulations. 
However there are few intervals with waiting time longer than r = 100^ d/g in K(t) because 
small noises chop up long intervals. Therefore applying the median filter of a time width of 
60 a/^/ ' g to K(t), we calculate the distribution of waiting time for the coarse-grained time 
series, p(t'). We find that p{r) decays approximately as a power function with D = — 2 as 
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Fig. 8. The distribution of the waiting time in K(t) for T = 2y / d/g and H = AOd. 
shown in Fig. 8. 

The plateau with a = —1.4 appears in a wide range of T as shown in Fig. 5(a), although 
the switchings can not be well-defined by K(t) as T increases. It is possible that the dynamics 
of the top location has some relations with the density fluctuation of flows. In experiments of 
granular flow in vertical pipes filled with fluid, the exponents —1, —4/3 and —3/2 are reported 
for the temporal power spectra of density, 22-24 ' 26-28,30 ' 31 and the exponents —1.4 and —3/2 
appear in traffic flows. 39-41 We note that these exponents are close to —1.4. 

As T becomes sufficiently large, we believe that a approaches —2. The top location stays 
at the same place for long time in comparison with T, and its displacements are caused by 
impulsive force of avalanches. We infer that the top location moves like as a random walk, 
and its power spectrum is Lorentzian-like, which decays as / -2 . 

5. Conclusions 

Carrying out 2-D DEM simulations, we have investigated the fluctuation of the top location 
of a sandpile that is caused by avalanches and piling up particles. We have found that the 
power spectra of the time series of the top location x top (t) behave as power functions in the 
range of long time scale. The exponent of the power spectrum, a, depends on the time interval 
T at which particles are fed to the sandpile. a is close to —1 for small T and decreases through 
a plateau with a = —1.4 as T increases. In the case of a = —1, avalanches occur mainly either 
on the left or right side slopes, and the states of the sandpile switch intermittently between the 
left and right modes. The power spectrum of the top location is approximately the same as that 
of the binarized time series defined from the switchings. In our simulations, the distribution 
of waiting time of the switchings obeys a power function with the exponent D = — 2 in this 
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case. The relation between D and a is consistent with the equation a + D = —3 proposed in 
the analytical theories. 37,38 
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